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Abstract 

Solving numerically hydro dynamical problems of incompressible fluids raises 
the question of handling first order derivatives (those of pressure) in a closed 
container and determining its boundary conditions. A way to avoid the first 
point is to derive a Poisson equation for pressure, although the problem of 
taking the right boundary conditions still remains. To remove this problem 
another formulation of the problem has been used consisting of projecting the 
master equations into the space of divergence free velocity fields, so pressure 
is eliminated from the equations. This technique raises the order of the differ- 
ential equations and additional boundary conditions may be required. High 
order derivatives are sometimes troublesome, specially in cylindrical coordi- 
nates due to the singularity at the origin, so for these problems a low order 
formulation is very convenient. We research several pressure boundary con- 
ditions for the primitive variables formulation of thermoconvective problems. 
In particular we study the Marangoni instability of an infinite fluid layer and 
we show that the numerical results with a Chebyshev collocation method are 
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highly correspondent to the exact ones. These ideas have been applied to lin- 
ear stability analysis of the Benard-Marangoni (BM) problem in cylindrical 
geometry and the results obtained have been very accurate. 



1 Introduction 

The problem of thermoconvective instabilities in fluid layers heated from below has 
become a classical subject in fluid mechanics |], ||. It is well known that two differ- 
ent effects are responsible for the onset of motion when the temperature difference 
becomes larger than a certain threshold: gravity and capillary forces. When both 
effects are taken into account the problem is called Benard-Marangoni (BM) convec- 
tion || . At the beginning theoretical studies considered layers of infinite horizontal 
extent without modelling lateral side-wall effects always present in experiments. 
More recently experiments have been concerned with confined containers with dif- 
ferent geometries [[§ 0, @, [7Jj an d consecuently theoretical and numerical studies 
have followed the same way [g, |], |T(J|, [□], [T2|, [T3], |T3 . 



Solving numerically hydrodynamical problems of incompressible fluids raises 
the question of handling first order derivatives for pressure in a closed container and 



determining its boundary conditions. In the velocity-pressure formulation [16|, [L7 



first order derivatives are avoided by deriving a Poisson equation for pressure [15 



although the problem of finding its boundary conditions still remains. Dirichlet 



and Neumann conditions are studied in detail in Refs. [f5, 17] and good results 
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are obtained with Neumann conditions in Ref. In this work, the sensitivity 

of the results, with respect to artificial pressure boundary conditions, is reported. 
The results we mention refer to rigid boundary conditions in the velocity field and 
purely hydrodynamical problems, but the issue of thermoconvection with free or 
Marangoni conditions for the velocity at the upper surface which we treat here has 
not been addressed. To avoid spurious results in the primitive variables formulation 
the finite difference and finite element methods use staggered grids, but this is not 
necessary with the collocation method we propose in this article. 

In thermoconvection problems pressure is usually avoided. In Refs. [B5, |20] 



the method of potentials of velocity is used to eliminate the variable of pressure 
from the equations. This technique raises the order of the differential equations and 
additional boundary conditions may be required. This is particularly troublesome in 
cylindrical coordinates where high order derivatives cause awkward difficulties. The 
spectral method used in Ref. [|14| allows the removal of pressure in the primitive 



variables formulation, however the collocation method that we study here is easier 
to implement. 

In this paper we present several pressure boundary conditions that allow us 
to solve thermoconvective problems in the primitive variables formulation. This is 
shown to be very useful when dealing with cylindrical coordinates. In the second 
section we perform the linear stability analysis of the Marangoni problem in an infi- 
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nite horizontal fluid layer. Since it has an exact solution it is used as a test problem 
to be compared with the numerical results obtained in the primitive variables for- 
mulation and several pressure boundary conditions. We use a Chebyshev collocation 
method and sensivity to the pressure boundary conditions is observed only in the 
convergence rate of the method. In the third section we study the BM problem 
in cylindrical geometry in the primitive variables formulation with the boundary 
conditions that provide better convergence for the test problem. We explain the 
methodology in detail. We study the convergence of the numerical results and com- 
pare them with the previous bibliography on the subject. In the fourth section 
conclusions are presented. 

2 Test problem 

The physical situation we consider corresponds to an infinite layer of a fluid with a 
free top surface. The domain is: 



In the vertical direction a temperature gradient is imposed. It is well known that 
motion sets in after the vertical temperature gradient has reached a critical value. 
The only mechanism that we consider in developing convection is the Marangoni 
effect, which consist of the variation of the surface tension with temperature. In 
the reference state there is no motion and heat propagates by conduction only. To 




(1) 



4 



study the linear stability of this state we must solve the linearized equations for 
the perturbations from the reference state. As there are no boundary conditions in 
the horizontal directions and the basic state is invariant under translations in those 
directions it is possible to take the Fourier modes in those coordinates. 



u' = u(z)e ikx + cc, 9' = 9(z)e ikx + cc, p' = p(z)e ikx + cc. (2) 

The linear problem to be solved is 

V-u = 0, (3) 
Pr _1 d t u = -Vp + Au, (4) 

d t 9 = u z + A9, (5) 

where u = (u x ,u y ,u z ) is the velocity field of the fluid, 9 is the temperature, p the 
pressure, V = (ik X) ik y , d z ), A = (d z — k 2 ), k = |k|, Pr is the Prandtl number which 
is considered infinite. The boundary conditions are: 

u|, =0 = 0, 9\ z=0 = 0, (6) 
(d z u x + Mik x 9)\ z=1 = 0, (d z u y + Mik y 9)\ z=1 = 0, (7) 
u z \ z=1 = 0, (d z 9 + B9)\ z=1 = 0, (8) 



where B is the Biot number and M is the Marangoni number which takes into 
account the variation of the surface tension with temperature. 
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2.1 Formulation without pressure 

Classically, pressure is eliminated and the following set of equations and boundary 
conditions are obtained, 

{D 2 -k 2 ) 2 u z = 0, (9) 

u z + {D 2 -k 2 )e = o, (io) 

u z \ z =o = Du z \ z=0 = 6\ z=0 = 0, (11) 
u z \ z=1 = (D 2 u z + k 2 M6)\ z=1 = (D6 + B6)\ z=1 = 0, (12) 

where D = d/dz. This formulation will be called PI in the following text. It can 
be easily shown that the u z and 9 solutions are 

u z (z) = ai e kz + a 2 ze kz + a 3 e~ kz + a A ze~ kz (13) 

, hy , hy 1 a\ze kz 1 ax z e kz 

0( z ) = b 1 e kz + b 2 e~ kz — 1 h 

v ; 2 k 2 k 

(lze kz lz 2 e kz \ (lze- kz lz 2 e~ kz \ 

° 2 4^~I^ + ° 4 4^ + I^ (14) 



Introducing them into the equations (p|)-(^2|) we get a linear system of equations for 
the unknowns a i; 6j. The solvability condition for the system gives the following 
dependence of the critical Marangoni number on k, 

8kB(a + Aak + a 2 - Aa 2 k - 1 - a 3 ) + 8k 2 (-l - a + Aka + a 2 + Aka 2 + a 3 ) 



M 



-3a 2 + Ak 3 a + 3a + Ak 3 a 2 - 1 + a 



where a is a function of k, a(k) = cosh(2 k) — sinh(2 k). 



2.2 Formulation in primitive variables 

We consider the primitive variables formulation with two different boundary condi- 
tions for pressure. 

1. Continuity equation at the boundaries 

(V-u)U=o,i = 0, (15) 

these boundary conditions together with equations (0)-© an d boundary con- 
ditions define the problem P2. 

2. Normal component of the Navier-Stokes equations at the top boundary and 
continuity equation at the bottom, 

(-d z p + Au z )\ z=1 = 0, (16) 
(V-u)U, = 0, (17) 

these boundary conditions together with equations (0)-© and boundary con- 
ditions ©-© define the problem P3. 

2.3 Numerical method 

The exact formulation (PI) together with formulations P2 and P3 are solved numer- 
ically with a Chebyshev collocation method. The eigenfunctions are approximated 
by Chebyshev polynomial expansions in the z direction. After changing the z coor- 
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dinate to transform the [0, 1] interval into [—1, 1] we use the expansions, 

N-1 N-1 N-1 

u = ]T a n T n (z), = £ b n T n (z), p=J2 CnTn(z), (18) 

n=0 n=0 n=0 

which are introduced into the equations. 

For problem PI the system and boundary conditions are evaluated at the 
collocation points, 

z i = cos (Ytt 1 ^ - 1 V) ( 19 ) 



k JV- 1 

In particular, (^) is evaluated at nodes i — 3, ...,N — 2; flTU|) at i — 2, ...,N — 1 



the first boundary conditions (|TT[) at i = N and the second ones (|12l) at i = 1. 
We obtain 2 x N unknowns and 2 x N equations. For problems P2 and P3: (|]) 
at the nodes i = 2, JV - 1; (|) at i = 2, JV - 1; © at z = 2, iV - 1; the 
first boundary conditions @ at z = iV; (^) and @ at i — 1. We obtain 5 x N 
unknowns and 5 x (JV — 2) + 8 equations. To complete the system we evaluate Eq. 

at i = 1, JV for problem P2, and Eq. @ at i = 1 and (0) at z = JV for 
problem P3. If the coefficients of the unknowns which form the matrices of the 
resulting linear system A and B satisfy det(v4 — XB) = 0, a nontrivial solution of 
the linear homogeneous system exists. This condition generates a dispersion relation 
A = X(k, M, B), equivalent to a direct calculation of the eigenvalues from the system 
AX = XBX, where X is the vector which contains the unknowns. When A becomes 
positive the basic state is unstable. In the critical situation A = with no imaginary 



part. In this case the dispersion relation can be written as M = M(k, B). We have 
calculated different critical Marangoni numbers for several k and B values. 

2.4 Convergence of the numerical method 

To carry out a test of the convergence of the method we compare the differences in 
the thresholds of M c for different orders of expansions. In table I the thresholds for 
k — 10 and B = 10 are shown for seven odd consecutive expansions. The thresh- 
olds converge to the exact value of 1600.01 and the difference between consecutive 
expansions tends to zero as N increases. Results on convergence improves greatly 
when k and B are smaller, i.e. for the critical value of k, k c as table II shows, for 
this reason there and in table III we use 7 polynomials expansions. 

2.5 Discussion 

As the exact solution is known, it is straightforward to compare the numerical meth- 
ods. In table II we show the critical Marangoni and wave numbers in each formula- 
tion for different Biot numbers. We calculate differences between each formulation 
and the exact solution. The mean relative error for PI and P2 is 10~ 2 , whereas the 
same for P3 is 2 • 10~ 3 . This indicates P3 as the best approach. In table III we 
show the critical Marangoni number for k = 10 in each formulation and different 
Biot numbers. Although convergence has not been reached (see table I) the result 
is significant because while the mean relative error for PI is 0.5 and for P2 it is 0.4, 
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for P3 it is 0.005. We confirm again that the third formulation converges better. 
Differences in accuracy between PI and P3 could be due to the fact that PI is a 
higher order problem. On the other hand P2 diminishes its efficiency with respect 
to P3 because of the boundary conditions considered. Although sensitivity of the 
numerical method to pressure boundary conditions has been addressed [jTj|, [L9j it 



seems here that it is not much more important than other effects such as increasing 
the order of the derivatives in the equations. Moreover we have shown that con- 
sidering appropriate boundary conditions for pressure leads to very good numerical 
performance. 

3 Benard-Marangoni in cylindrical geometry 
3.1 Formulation of the problem 

The physical situation corresponds to a fluid layer of thickness d filling a cylinder 
with radius I. The surface tension at the upper free surface is temperature dependent 
and the fluid is heated from below. Motion sets in when the vertical temperature 
gradient has exceeded a critical value. 

In the basic state, there is no motion and heat propagates by conduction. In 
the linear Boussinesq equations, the field can be expanded in the azimuthal variable 
4> in Fourier modes as follows, 

u(t,r,<f>,z) = u m (t,r,z)e m ^ (20) 
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6(t,r,<f>,z) = 9 m (t,r,z)e im *, (21) 
p(t,r,cf>,z) = p m (t,r,z)e im ^ (22) 

where 9 and u = ue r + ve^ + we z are the infinitesimal temperature and velocity 
perturbations with respect to the conductive solution, p is the pressure perturbation, 
and (r, 0, z) are the polar coordinates. If space, time, velocity, temperature and 
pressure fields are respectively divided by the constants d, d 2 /n, n/d 2 , AT and 
PoKv/d 2 , the equations in dimensionless form are obtained. Here k is the thermal 
diffusivity, po is the mean density, rj is the dynamic viscosity (related to the kinematic 
viscosity through the expression, v = rj/po) and AT is the (conductive) temperature 
drop between the bottom and the top layers. The linearized general equations for 
the perturbations of index m are easily shown to be 



„ ! dv im - 2im , n4 . 

Pr ft = -— P + V ^ + lV + —U, (24) 

Pr-^ = -f z+ V 2 m2 w + R9, (25) 
f)9 

g~ t = w + V 2 m ,9, (26) 

= 1^)+^ + ^, (27) 
r or r oz 

where the subindex m has been cancelled for simplicity. The Rayleigh number is 
defined by R = agATd 3 / up, where a is the coefficient of volume expansion and 
g is the gravity. The Prandtl number is given by Pr = v/k. is defined by 
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V 2 n = r~ l d/dr(rd/dr) — nr~ 2 + d 2 /dz 2 . 

The boundary conditions are the following. The bottom of the box is rigid 

and assumed to be perfectly heat conducting so that 

~ „ ~ 1 d(ru) im dw n „ 
u = 0, = 0, — + — v + tt- = on 2 = °- (28) 
r or r oz 

The upper surface of the fluid is assumed to be plane, nondeformable and free where 
surface tension effects are taken into account. We also assume that at the top, heat is 
transferred from the liquid to the ambient gas according to Newton's law of cooling, 
which results in a Biot condition for the temperature perturbations. For pressure 
the normal projection of the Navier-Stokes equations into the top plane is considered 
as a boundary condition. The mathematical expressions of the boundary conditions 
at the upper surface are then, 

w = 0, ^ + B6 = 0, ^ + M^ = 0, ^ + M^ = 0, on z = 1, (29) 
oz oz or oz r 

-^- + V 2 m2 w + R6 = 0, on z = 1, (30) 
oz 

where B is the Biot number and M = jATd/ PqKU is the Marangoni number with 
7 the constant rate of change of surface tension with temperature. 

The lateral side wall is rigid and we will consider it adiabatically insulated. 
For pressure we consider the normal projection of the Navier-Stokes equations into 
this wall as boundary condition. This is written as follows 

u = 0, 86/ dr = 0, on r = a, (31) 
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dj) 27/77 

where a = l/d is the aspect ratio. 

The use of cylindrical coordinates, which are singular at r = 0, imposes 



regularity conditions on the velocity and temperature fields ||22|| . These conditions, 
which express that the unknown fields are single valued at r = 0, are mathematically 
summarized as follows, 

where x is any scalar or vectorial field. For the temperature and velocity perturba- 
tions, these conditions are written as 

dw 88 dp , . 

n= 7r = 7T = 7r = ' m = ' 34 
or or or 

u + iv = w = 9 = p = 0, m=l, (35) 
u = v = w = 9 = p = 0, m^O, 1. (36) 
3.2 Numerical method 

The numerical method used is a collocation method which is similar to the method 



used in Ref. [23| for rectangular containers. Firstly we change the variables z and 
r to transform the intervals [0,1] x [0, a] into [—1,1] x [—1,1] to use the Cheby- 
shev polynomials. The fields in equations (CT)-(B7p, denoted generically by x, are 
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aproximated by Chebyshev expansions, 

L-l AT-1 

* = EE< T 'M T »W. (37) 

Z=0 n=0 

which are introduced into the equations. The resulting expressions are evaluated at 
the collocation points (rj, Zi) defined as follows, 

r 3 = cos (Q-^i - l) tt) , (38) 

* = cos ((^T~TI ~ 1 ) 71 ) ' (39) 

In the case m > 1 the system and boundary conditions are evaluated at the following 
collocation points: Eqs. (|23| ) -(|27| ) at the nodes i = 2,...,N — 1, j = 2,...,L — 1; 



the boundary conditions at z — — 1 ( |28|) at i — 1, j = 2, ...,L — 1; the boundary 
conditions at z = 1 © at i = JV, j = 2, L — 1 and ® at i = N, j = 2,..., L; 
the boundaries at r = — 1 ( P^) at z = 1, iV, j = 1; finally the boundaries at r = 1 
@ at z = 1, JV, j = L and (H) at z = 1, JV - 1, j = L. We obtain 5xiVxI 
equations and 5 x N x L unknowns. For m — 1 the evaluation has been done at 
the same nodes, but as there are only four boundary conditions at r = — 1 we have 
diminished the expansion range of the v field by one order, 

L-2N-1 

, = EEi T 'M T «(4 (4°) 

1=0 n=0 

Therefore we get 4xNxL + Nx (L — 1) equations with the same amount of 
unknowns. In the case m = the angular component of the velocity is nulle and Eq. 

14 



(|24|), together with the corresponding boundary condition for v in Eqs. (p8|)-(p6|) 
disappear and therefore we obtain 4 x N x L equations with the same amount of 
unknowns. In this case the matrix associated to the linear algebraic system is sin- 
gular, due to the fact that pressure is defined up to an additive constant. To fix this 
constant in the node i = N — 2, j = L the boundary condition (|32|) is replaced by a 
Dirichlet condition for pressure (i.e, p = 0ati = N — 2,j = L). The resulting linear 
systems have been solved with a standard numerical package. The eigenfunctions 
and thresholds of the generalized problem are numerically calculated. If the coeffi- 
cients of the unknowns which form the matrices A and C satisfy det(A — AC) = 0, 
a nontrivial solution of the linear homogeneous system exists. This condition gen- 
erates a dispersion relation A = X(R, M, B), equivalent to a direct calculation of 
the eigenvalues from the system AX = XCX, where X is the vector which contains 
the unknowns. When A becomes positive the basic state is unstable. In the critical 
situation Re(X) = without imaginary part. In this case the dispersion relation can 
be written as M = M(R, B) or R = R(M, B) and critical Marangoni number or 
critical Rayleigh number are obtained directly from the eigenvalue problem. 

3.3 Convergence of the numerical method 

To carry out a test of the convergence of the method we compare the differences in 
the critical Marangoni number M c to different orders of expansions for different val- 
ues of the aspect ratio a. In table IV these thresholds are shown for four consecutive 
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expansions varying the number of polynomials taken in the z direction (N) and in 
the r direction (L). These results allow us to conclude that N x L = 9 x 13 gives 
very good results for aspect ratios lower than 5. So, except where otherwise stated, 
all results given below correspond to the values N x L = 9 x 13. 

3.4 Comparison with other theoretical works 

Our main goal is to show how a Chebyshev collocation method applied to the prim- 
itive variables formulation of thermoconvective problems with convenient boundary 
conditions for pressure produces excellent results. For this reason we compare our 
results with the preceding ones on the subject. 

Dauby et al. have already considered the comparison between their results 
and the preceding ones. For this reason we are focusing on a comparison only with 
their paper. The correspondence between our work and that of Dauby et al. can 
be tested by comparing our Tables V and VI to their Tables II and III. ^From this 
comparison we see that the deviations are less than 0.06%. We have checked that the 
succession of unstable eigenmodes when a is increasing in our approach coincides 
with those of Dauby et al.. In their work they discuss differences obtained with 
Zaman and Narayanan ||13|| . As is shown in figure 1 when a increases all the critical 
values converge to the same thresholds and it is harder to distinguish among them 
with a numerical method. 
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4 Conclusions 

We have shown that a Chebyshev collocation method applied to the primitive vari- 
ables formulation of thermoconvective problems with convenient boundary condi- 
tions for pressure produces excellent results. First, we have solved the Marangoni 
problem in infinite geometry with a Chebyshev collocation method, keeping the 
primitive variables formulation and testing two sets of boundary conditions for pres- 
sure. We have compared these results with the exact ones obtained classically by 
eliminating pressure. As there is an exact solution to this problem it was straight- 
forward to compare the performance of the numerical results in the different for- 
mulations. We prove that the results obtained do not depend on those boundary 
conditions. Sensitivity to those boundary conditions is only observed in the con- 
vergence of the method. The best results are obtained in the primitive variables 
problem, taking as boundary condition for pressure the projection of Navier-Stokes 
equations at the top boundary and continuity at the bottom one. Second, consid- 
ering the boundary condition with the best convergence we have solved the BM 
problem in cylindrical geometry. We have compared with previous numerical re- 
sults and the correspondence is very high. We conclude that the primitive equations 
with the appropriate boundary conditions for pressure of these thermoconvective 
problems are accurately solved with our collocation method. 
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Table I 



N 


PI 


P2 


P3 


17 


1600.01 


1600.01 


1600.01 


15 


1599.99 


1599.91 


1600.00 


13 


1599.42 


1597.95 


1599.66 


11 


1590.73 


1572.61 


1595.49 


9 


1500.53 


1382.53 


1572.57 


7 


972.58 


776.62 


1622.50 


5 


231.32 


201.84 


2817.82 



Table II 



B PI (numerical) 


PI ( 


exact) 


P2 


P3 


M c k c 


M c 


k c 


M c k c 


M c k c 



0.1 83.20 2.04 83.43 2.03 82.93 2.05 83.31 2.03 
1 115.69 2.26 116.13 2.25 115.14 2.28 115.92 2.25 
10 410.45 2.77 413.44 2.74 406.27 2.82 412.20 2.75 



Table III 



B 


PI (numerical) 


PI (exact) 


P2 


P3 


0.1 


484.39 


808.01 


386.80 


808.08 


1 


528.77 


880.01 


422.23 


882.12 


10 


972.58 


1600.01 


776.62 


1622.50 



Table IV 



m 




5x9 


7 x 11 


9 x 13 


11 x 15 


m = 





150.705 


153.689 


154.064 


154.062 


m = 


1 


148.156 


152.558 


152.945 


152.949 


m = 


2 


150.946 


153.734 


154.154 


154.154 


m = 


3 


148.526 


152.899 


153.255 


153.256 


m = 


4 


149.400 


153.542 


153.979 


153.987 



Table V 
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a = 


1 


a = 


I 


a 


= 4 






M c 


M c 


M c 


M c 


M c 


M c 


B(R = 


o) 


(H. and M.) 


(D. et al.) 


(H. and M.) 


(D. et al.) 


(H. and M. 


) (D. et al.) 


0.01 




164.65 


164.55 


84.640 


84.638 


82.485 


82.486 


0.1 




168.44 


168.34 


88.346 


88.344 


86.263 


86.263 


1 




206.29 


206.16 


125.011 


125.000 


120.624 


120.630 






R c 


R c 


R c 


R c 


R c 


R c 


B(M = 


0) 


(H. and M.) 


(D. et al.) 


(H. and M.) 


(D. et al.) 


(H. and M. 


) (D. et al.) 


0.01 




1419.30 


1419.47 


712.542 


712.667 


695.543 


695.668 


0.1 




1426.06 


1426.24 


726.574 


726.704 


709.326 


709.457 


1 




1481.89 


1482.12 


835.897 


836.072 


799.522 


799.735 



Table VI 



a = 1 a = 2 a = 4 a = 8 

m = 163.676 80.878 78.777 76.179 

m = 1 108.383 91.254 77.864 76.448 

m = 2 158.994 98.407 79.699 76.203 

m = 3 255.885 99.955 78.095 76.487 
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TABLE CAPTIONS 
Table I 

Critical Marangoni number M c for different orders of expansions and B = 10, k = 
10. 

Table II 

Critical Marangoni number M c and wave number k c for different values of the Biot 
number B. 

Table III 

Critical Marangoni number M c for different values of the Biot number B and k = 10. 
Table IV 

Critical Marangoni number M c for different orders of expansions and B = 2, a = 5. 
Table V 



Comparison between our results (H. and M.) and those of Dauby et al. [ I4| . The 
first lines give the critical Marangoni number when R = and for different Biot 
numbers; the bottom of the table presents the critical Rayleigh number for M — 0. 

Table VI 

Critical Marangoni numbers M c for different azimuthal wave number m. The 
Rayleigh and Biot numbers are 100 and 0.2, respectively. 
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FIGURE CAPTIONS 
Figure 1 

Critical Marangoni number M c as a function of the aspect ratio a. Curves corre- 
sponding to different azimuthal wave numbers m are represented. The Rayleigh and 
Biot numbers are zero. 



25 



